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Abstract 

Several algorithms for introducing artificial dissipation into a central difference approx- 
imation to the Euler and Navier Stokes equations are considered. The focus of the paper is 
on the convective upwind and split pressure (CUSP) scheme, which is designed to support 
single interior point discrete shock waves. This scheme is analyzed and compared in detail 
with scalar and matrix dissipation ( MATD ) schemes. Resolution capability is determined 
by solving subsonic, transonic, and hypersonic flow problems. A finite-volume discretiza- 
tion and a multistage time-stepping scheme with multigrid are used to compute solutions 
to the flow equations. Numerical results are also compared with either theoretical solutions 
or experimental data. For transonic airfoil flows the best accuracy on coarse meshes for 
aerodynamic coefficients is obtained with a simple MATD scheme. 


*This research was supported in part by the National Aeronautics and Space Administration under NASA 
Contract No. NASI- 19480 while the first author was in residence at the Institute for Computer Applications in 
Science and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23681-0001. 
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Introduction 


Accuracy must be a primary consideration in the construction of any numerical scheme. One 
would like to devise a scheme with the minimum amount of artificial dissipation required 
for stability, as well as convergence in the case of a stationary solution. For fluid dynamic 
computations the numerical scheme should be designed to have high accuracy in smooth regions 
of the flow field and well resolved shock waves and contact discontinuities. According to 
Harten [3] such discrete formulations, where the accuracy away from discontinuities is at least 
second order, are called high resolution schemes. The design of these schemes for systems of 
conservation laws is generally based on theory developed for a scalar conservation law. As a 
consequence one cannot ensure that the properties of the scheme for the scalar equation are 
valid for the system. In addition, schemes that permit high definition of shock waves without 
oscillations are first order in the neighborhood of shocks. Concern naturally arises regarding 
contamination of the solution, especially in the case of viscous flows. For these reasons the 
properties and resolution capability of this class of schemes must be verified through numerical 
applications for a wide range of flow conditions. 

High resolution schemes of particular interest for solving the compressible fluid equations 
are those that allow shock capturing with a single interior point. In Ref. [ 6 ] Jameson presents 
the convective upwind and split pressure (CUSP) scheme. For this scheme the artificial diffusive 
flux vector associated with a given coordinate direction is expressed in terms of changes in the 
state and flux vectors. A somewhat limited number of inviscid and viscous computations have 
been performed to evaluate these schemes (see Refs. [ 6 , 7] and [22, 23, 9]). 

In the current work we investigate and analyze the HCUSP version [ 6 ] of the CUSP scheme. 
The HCUSP scheme allows a solution with constant total enthalpy for steady flow. We discuss 
the shock-capturing behavior for various choices of the dissipation coefficients. We introduce 
a simple modification of the limiter function, which is generally used with the scheme, to 
control background dissipation, and thus global accuracy. The HCUSP scheme includes a 
contribution that is scaled according to the local velocity. If the velocity approaches zero, 
as it does for viscous flow near a solid boundary, and there is a high aspect ratio mesh, 
the dissipation in the streamwise direction (i.e., direction of long side of mesh cell) may not 
be adequate for convergence. A change in the velocity scaling factor based on aspect ratio is 
presented. The resolution capability of the HCUSP scheme is evaluated for subsonic, transonic, 
and hypersonic flow problems. A detailed comparison of the scheme with a scalar and matrix 
dissipation schemes is performed. The scalar scheme is based on the dissipation model of 
Jameson, Schmidt, and Turkel [4] and is often used with central differencing. 


Dissipation 

A finite- volume approach is applied to discretize the fluid dynamic equations of motion. The 
computational domain is divided into quadrilateral cells, fixed in time, and for each cell the 
governing equations can be nondimensionalized and written in integral form as follows: 

~ [[ wdxdy + t ( fdy - gdx) = f ( f v dy - g v dx), ( 1 ) 

at JJu Jdn He Jdu 

where Q is a generic cell with boundary Sfi. In the scaling factor for the viscous terms on the 
right hand side of Eq. ( 1 ), the quantities 7 , M, and Re are the ratio of specific heats, Mach 
number, and Reynolds number, respectively, with M and Re defined in terms of nominal con- 
ditions. Taking wj ^ as the cell-averaged solution vector, Eq. (1) can be written in semidiscrete 
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form as 

+ £ w j,k — o, 

where Ctj } k is the area of the cell, and C is defined by 

C = £c + C D + Cad, 

with the subscripts C, D, and AD referring to convection, diffusion, and artificial dissipation. 
In order to simplify the description of the dissipation models, we consider the one-dimensional 
Euler equations of gas dynamics. 


Scalar Dissipation Model 

The scalar dissipation is based on the model introduced by Jameson, Schmidt, and Turkel 
[4]. This model defines a switching function based on a blending of the second and fourth 
differences. The term associated with the operator Cad is expressed as 


Then 


Cad Wj = ~{D 2 - D i )wj = d j+1/2 - d 3 _ 1/2 . 


D 2 Wj 

D 4 Wj 


V 

V 


(\'+ 1/2 

(-S+1/2 


£ \+l/2)^\ W 3' 

$1/ 2 )AVA] wj , 


( 2 ) 

( 3 ) 


where the index j refers to a cell center, and the operators A and V are forward and backward 
difference operators. The variable scaling factor A j+i/ 2 is defined as 


Aj+1/2 - 2 + (\)j+iL 


where A^ is the largest eigenvalue in absolute value (i.e., spectral radius) of the flux Jaco- 
bian matrix associated with the the Euler equations. In two dimensions, with (£, 77 ) denoting 
arbitrary curvilinear coordinates, the scaling factor is usually defined as 

= < t ) j,k{ r ) (A = 1 + r ^k 

with r = The exponent £ is generally taken to be between \ and |. The coefficients 

and £ use the pressure as a sensor for sharp gradients, and they are defined as 


£ 5 +i/2 = * (2) vj, v j+u v j+2 ), 

v _ Pj-i ~ 2 Pj + Pj+i 

3 pj-i + 2pj + p J+l 

4+i/a = majc [°> ( k<4) _ ^+ 1 / 2 )] • 


( 4 ) 


where typical values for the constants and are in the ranges j to ~ and ^ to 
respectively. We shall refer to Eq. (4) as the JST switch. The switching function v can be 
interpreted as a limiter, in the sense that it maximizes the second-difference contribution at 
extrema and switches off the fourth-difference term. Moreover, at shock waves the dissipation 
is first order, and a first-order upwind scheme is produced for a scalar equation. In smooth 
regions of the flow field the dissipation is third order. 
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Thus, we have two different dissipation mechanisms at work. The switch determines which 
one is active in any given region. For smooth flows, v is small and the dissipation terms consist 
of a linear fourth difference that damps the high frequencies which the central difference scheme 
does not damp. This is necessary for achieving a steady state. In the neighborhood of large 
gradients in the pressure, v becomes large and switches on the second-difference viscosity 
while simultaneously reducing the fourth-difference dissipation. This is needed to introduce 
an entropy condition so that the correct shock relations are satisfied and to reduce overshoots 
near discontinuities. 


Matrix Dissipation (MATD) 

Sharp resolution of shock waves without oscillations can be achieved by closely imitating an 
upwind scheme in the neighborhood of a shock wave. A key feature of upwind schemes is a 
matrix evaluation of the numerical dissipation. With this evaluation the dissipative terms of 
each discrete equation are scaled by the appropriate eigenvalues of the flux Jacobian matrix 
rather than by the spectral radius, as in the JST scheme. A matrix dissipation model can 
easily be constructed by starting with the JST formulation. 

One can show [20] that the necessary modification of the JST scheme to produce a matrix 
dissipation model is the substitution of \A\ for the eigenvalue scaling factor A in Eqs. (2) and 
(3). Since the Euler equations are a strongly hyperbolic system, the coefficient matrix can be 
diagonalized. Assume QAQ~ l = A (diagonal matrix). Then \A\ is defined as \A\ = Q~ l \A\Q 
and |A| = diag{\\\ |, | A2 1 , | A3 1 ) , where \ are the forward acoustic, backward acoustic, and 
convective eigenvalues. Efficient ways of computing |A| times a vector axe presented in Ref. [20]. 

In practice one cannot choose Ai,A 2 ,A 3 as the eigenvalues. Near stagnation points A3 
approaches zero while near sonic lines Ai or A2 approaches zero. A zero artificial viscosity 
would create numerical difficulties. Hence, we limit these values as 

A* I = max(|A i |,V r n />(yl)), j = 1,2 
M = max (|A 3 |, V t p{A)), 

where p(A) is the spectral radius of A, and the linear eigenvalue A3 can be limited differently 
than the nonlinear eigenvalues. The parameters V n and Vi have been determined by numerical 
experimentation. Typical values axe V n = 0.25 and Vi = 0.025. 


TVD and LED Properties 


In one dimension consider the approximation 

w j +1 = w i ~ 2 Ax ^ j+1 ~ + 2 Ax ^ w j + i /2 - Qj-i/2& w j-i/2)i 

where A u?j +1 /2 — Wj + 1 — Wj- Suppose that the Jacobian matrix A = is at least locally 
constant. Then the scheme is TVD if 


Qj+ 1/2 > 


A j+ 1/2 ■ 


(5) 


This result follows directly from the fact that the system can be diagonalized, and we can 
ensure that each characteristic variable satisfies the TVD criterion (i.e., the total variation of 
the variable is nonincreasing). 
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A replacement for the switch of Eq. (4) that is TVD for a scalar equation is introduced in 
Ref. [20]. In one dimension this switch is given by 


_ jgj+i - +pj- 1| 

3 \Pj+i ~Pj\ + \Pj -Pi-l | + e' 


(6) 


and choose k ^ 


In practice we usually use a weaker form than Eq. (6), for example, 


_ \Pj+i ~2pj +Pj-i[ 

J (1 — u)Vtvd + wP ’ 


where 


■PrvD = |Pj+i - Pj| + | Pj ~ Pj - i |, 

V = p j+ i + 2 pj + pj. i, 

and 0 < w < 1. The TVD switch of Eq. (6) is recovered when Typically a; 1/2. 

In Refs. [6] and [8] Jameson develops a theory for scalar nonoscillatory schemes based on 
the local extremum diminishing (LED) principle that maxima do not increase and minima do 
not decrease. The LED principle applies to multidimensional problems, and it is equivalent to 
the TVD principle in one dimension. If a scheme is LED , then the scheme is positive. 


SLIP Scheme 


Until now we have considered a combination of a low-order and high-order artificial viscosity 
based on a scalar switch. This switch has the disadvantage that it is based on only one quantity, 
the pressure. Moreover, it forces all variables to be treated equal, even though some experience 
sharp changes through the discontinuity while others are continuous across the shock. One 
can instead limit independently each dependent variable in each coordinate direction. 

Jameson [5] introduced a new class of limiters and implemented such a limiting process 
within a framework that he called the symmetric limited positive (SLIP) scheme. For the SLIP 
formulation Jameson constructed a family of limiter functions based on the function 


R(u : v) = 1 — 


u — v 


<? 


u\ + |v| + e 


where q is a positive number and e has the dimensions of u. Note that i?(u, t>) « 0 whenever u 
and v have the opposite sign. Let w be an element of the solution vector for the flow equations. 
According to our previous theory [20] R(Awj +3 / 2 ^ w j - 1 / 2 )^ where AtVj+ 3 / 2 = w j + 2 — Wj+ 1 , 
would be replaced by where ^*+ 1/2 is the maximum of Vj over the nearest neighbors 

and v is given by Eq. (6). Define the limiter function L(u, v) by 

L(u, v ) = R(u, (7) 

At the mesh cell interface j + 1/2, we define the left and right states for each dependent variable 
as 

WL = Wj + ):L{ AtUj+3/2, A Wj_ 1/2 ), (8) 

w R = w j+ i - ^L{Aw j+3/2 ,Awj_ 1/2 ), 
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and so 


U)R-W L = Aw j+l/2 - L(Aw j+3/2 , Awj-ip). 

For the artificial viscosity all differences will be based on wr — wl- In the neighborhood of 
shock waves R(u,v) and hence are close to zero. Moreover, wr — wr — 

resulting in a first-order scheme for the artificial viscosity. For smooth flow R(u,v) ~ 1, and 
L(u, v) = (u + v)/2. Hence, in a smooth region 


W R -W L = Awj +1 / 2 - L(Aw j+3 / 2 , Awj_ l/2 ) 

A Aw j+ 3/2 + AWj _\/2 

- &Wj+ 1/2 2 ( 9 ) 

1 a 3 

= - -A W j+ 1/2 . 

Thus, the SLIP scheme has similar properties to the JST scheme. One can obtain the rela- 
tionship between the SLIP and JST schemes by defining the diffusive flux dj+ 1/2 as 


dj+ 1/2 = a j+1/2 (w R - w L ) (10) 

with ctj+ 1/2 = K' 2 ® A j+1 / 2 . The quantity is a parameter, and A is the spectral radius of the 
associated flux Jacobian matrix. 

One difference between the JST and SLIP schemes involves the parameters k ^ and 
for the second and fourth differences, respectively. Both and are free parameters in 
the JST scheme. As seen from Eqs. (10) and (9) these parameters are automatically chosen as 
and with the SLIP scheme. The coefficient of the second difference is chosen as \ 

so that the scheme is fully upwind for supersonic flows. However, the purpose of the fourth- 
difference viscosity is only to accelerate the convergence to a steady state by eliminating the 
odd-even point decoupling. Hence, we wish to be as small els possible for accuracy while still 
achieving a good convergence rate. It does not seem reasonable to connect the two components 
of the artificial viscosity. 

One can generalize the SLIP scheme by reintroducing a free parameter that essentially 
governs the level of the third-order viscosity in smooth regions. The resulting scheme has the 
disadvantage that a free parameter must be chosen; however, it has the advantages of greater 
flexibility and increased accuracy. We replace Eq. (7) by 


L(u, v,w) = i?(u, w) * 


(1-4k(> + W 4 >^1, 


( 11 ) 


where left and right state values arc determined by 

w L = Wj + ^L(Aw j+3/2 , Aw j+1/2 , Awj_i/ 2 ), 

WR = w j+ 1 - ^-L(Aw j+3/2 , Aw j+1/ 2, Awj_ l/2 ). 

When k) 4 - = 1, the L of Eq. (11) reduces to the original L of Eq. (7). Now, at shock waves 
R(u,w) is small and we have wr — wr = A Wj+ 1 / 2 * For smooth regions of the flow field we 
have wr — wl = — 2k^ A 3 u^ +1 / 2 - 
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CUSP/AUSM Scheme 

So far we have described the use of an artificial viscosity based on either a scalar or ma trix co- 
efficient. Liou and coworkers designed a scheme called Advection Upstream Splitting Method 
(AUSM) [10, 11, 12]. This method was later refined for large-scale 3-D viscous computations 
[15]. AUSM is based on a splitting of the flux function into convective and pressure contribu- 
tions. In some sense, the pressure terms contribute to the acoustic waves while the velocity 
terms contribute to convective waves. Hence, it is reasonable that these flux terms be treated 
differently. Liou thus considers decompositions of the flux vector that are not based on a char- 
acteristic decomposition but on Mach number scaled contributions of the left and the right 
states to the interface flux. This decomposition has the disadvantage that it is more difficult 
to develop for other sets of equations compared with a characteristic decomposition. A similar 
type scheme called the convective upwind split pressure (CUSP) scheme was later introduced 
by Jameson [5] and subsequently modified by Tatsumi, Martinelli, and Jameson [7, 8, 22, 23]. 
The CUSP scheme has several advantages. First, one can consider the scheme as another type 
of artificial viscosity, since it is defined as a sum of the central flux average plus a dissipative 
flux. Hence, it can be readily used with a variety of time-stepping schemes (e.g., multistage, 
LU, implicit, etc.). Second, the CUSP formulation can be used with multistage schemes which 
do not evaluate the artificial dissipation fluxes at every stage, in order to reduce computa- 
tional work. Another advantage of the CUSP scheme is that it can be easily combined with 
preconditioning, since preconditioning is based on the inviscid flux form and not the artificial 
dissipation. Hence, we shall only describe the CUSP version of this type of scheme. 


Definition of CUSP Scheme 


Previously, we introduced the scalar and matrix- valued viscosities by considering cL+ 1/2 of the 
form 

d j+l/2 = 2<?i+l/2(%l - Wj ). 

The factor \ is introduced so that we get first-order upwinding when Qj+i /2 = I ■ We note 
that for the scheme to be TVD, Q must be sufficiently large (see Eq. (5)). For the matrix 
viscosity we chose Q = |A| (modified near zero eigenvalues) while for the scalar viscosity we 
chose Q — cr(A)I. 

For the CUSP scheme we instead choose d as a linear combination of w and /. We shall 
only consider the choice for the state vector given by 


w h = (p pu pH) 7 


( P 

f — u pu 

\ pH 


+ 


This choice is denoted HCUSP by Jameson 
defined as 


( ° 

P 

v° 

[7]- 


j = uw h + / p . 


The first-order accurate CUSP scheme is 


d i+ 1/2 = - Wj) + ^(fj+1 - fj). (12) 

The factor c is included so that v is dimensionless. We thus consider only scalar parameters 
instead of a matrix coefficient, but we have two free parameters, v and /?. The scheme is total 
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enthalpy preserving. By using the arithmetic average, u — ^(uj+i + Uj), and the definition 


otc— vc 4- (3u , 


( 13 ) 


one can rearrange Eq. (12) to obtain 

d j+ 1/2 = ^otc(w j+ 1 - Wj) + ^(fp j+1 - f Pj ) + ^w( u j+i - u j )• 

Introducing the Roe matrix [17] Arl , we have that Jr — Jl = A/^iur — ^l)- This relation 
is exact if Arl is computed from weighted averages of the left and the right states. Then the 
first-order dissipation is 


dj+1/2 = ^(P A RL + vcI)(w R - W L ). 


(14) 


We see from this formula that d is a linear function of A. Recall that | A| is a quadratic function 
of A, by the Cayley-Hamilton Theorem. Hence, it is not possible to bound d by \A\. Thus, 
the CUSP scheme cannot be TVD. 

Assume that the subscript L denotes the interior point inside the shock zone, R is the 
state downstream of the shock, and the state LR is subsonic. Jameson [7] shows that the 
downstream point with the state R is in equilibrium if 


fn~ fL+ 1 + p ( w R ~ W'l) = 0, 


(15) 


Substituting the Roe matrix for the difference in / into Eq. (15) we get 


^Arl + •-I'j ( wr - wl) = o. (16) 

Hence, wr — wl is an eigenvector of Arl , and — (z^c)/(l + (3) is the corresponding eigenvalue. 
The eigenvalues of Arl are known to be A + , A~ and u. If A is an eigenvalue of A, then using 
this formula for vc in Eq. (14) we have 

d j+ 1/2 “ 2 [-XI + (3 {A rl - XI)] {w R - w L ). 

In order to have a positive diffusion when u > 0, we require that A be negative (i.e., — (i/c)/(l + 
(3) — A - ). Thus, 


vc = —{l + (3)X . 


(17) 


For u < 0 we obtain similarly 


vc= {1 — /?)A + . (18) 

So, we have reduced our two free parameters to one free parameter by demanding a one point 
shock profile. More generally, Jameson shows that one obtains a shock profile with one interior 
point if the following two conditions hold: 
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I. When the flow is supersonic through the shock then one obtains a totally upwind flux. 
II* The artificial dissipation Q satisfies a generalized eigenvalue problem 


(Arl ~ chraQra) (®a — w A ) = 0 

at the exit from the shock. 


The second condition is satisfied by both the matrix viscosity and the CUSP scheme^ however, 
the scalar viscosity does not satisfy the first condition. We again note that the TVD condition 
Q > \A\ is satisfied by the scalar and matrix viscosities but not by CUSP. 

What remains to be done is to choose suitable functions for (3 and vc which satisfy the 
above requirements. Jameson s choice for /3, which is based upon the eigenvalues corresponding 
to the acoustic waves, is given by 



+ maxfo,“+M 

if 

0 < M < 1 


(3=< 

- max ( 0 , j 

if 

-1 < M < 0 

(19) 


, sgn(M) 

if 

\M\ > 1. 



The cutoffs, (3 > 0 for u > 0 and P < 0 for u < (). ensure that the pressure terms are discretized 
centrally for small Mach numbers. Shock capturing with one interior point is obtained by taking 


vc = < 


-(1 + /?)A- 
(1-/?)A+ 
0 


if (3 = 0 

if (3 > 0, 0 < M < 1 
if (3 < 0, -KM <0 
if \M\ > 1. 


The dissipation coefficients are to be computed from Roe-averaged quantities. They provide 
full upwinding for supersonic flow {(3 = sgn(M), i/ = 0). The choice of vc = |u| for /3 = 0 yields 
a continuous dissipation coefficient in the subsonic region, and it does not smear slip lines with 
|u| close to zero. This makes the CUSP formulation attractive for viscous flow calculations with 
boundary layers. However, viscous flows are usually discretized by using cells with large aspect 
ratio. It is well known that this situation requires larger dissipation scaling in the direction of 
the long cell sides than given by |u|. We redefine the dissipation coefficients in the individual 
coordinate directions. For the f -direction we have 


( max(|u|, 6cr ) if 
— (1 + (3) A if 


vc^ = r < 


(1 - (3) A+ if 

0 if 


(3 = 0 

(3 > 0 and 
0 < M < 1 
(3 < 0 and 
-1 < M < 0 
\M\ > 1. 


(20) 


where r + and r are functions of the spectral radii in the £ and rj directions (A ? and A^), and 
they are defined as follows: 


r = (T^) C , r + = max(r, 1), 
a € 


r = min(r, 1). 


The dissipation coefficient in the 77 -direction is defined correspondingly. 


8 



Simplified Scheme 

Several modifications of the CUSP scheme have been in use so far. Based upon the w ^ = 
(p pu pH) T system the dissipation coefficients presented in Refs. [7] and [23] are as follows: 


1 

r 

\M\ if 

m 

> 6 

a = < 

1 1 

i‘+*r) « 

\M\ 

< e, 


( 

+ 


) if 

0 < M < 1 

(3 = < 

- 

max (0, 

) if 

-1 < M < 0 


k. 

sgn(M) 

if 

\M\ > 1. 


This choice does not allow exact shock capturing because Eqs. (17) and (18) are not satisfied. 
Furthermore, Roe averaging has been replaced by arithmetic averaging in Ref. [7] and A”, A + 
by u — c, u + c, respectively. This simplification saves a few square roots in the coding of the 
dissipative flux. This then becomes 


a 


\M\ if \M\ > e 
1( C+ M!) if |M| < e, 


0 


max (0, 2 M — 1) 
< min (0, 2M + 1) 
sgn(M) 


if 0 < M < 1 
if -1 < M < 0 
if \M\ > 1. 


( 22 ) 


Higher Order Scheme 

Having determined vc and / 3 , we see from Eq. (12) that the scheme is completely defined in 
terms of w and /. Recalling from previous sections, we now need to combine a first-order 
artificial dissipation with a high-order dissipation. As previously the high-order dissipation is 
used in the smooth regions w T hile the low-order artificial viscosity is used in the neighborhood 
of shocks. To allow changing from one type to the other we can either use the JST switch of 
Eq. (4) or the SLIP scheme. Formula (12), as given, is only first-order accurate, as it depends 
only on dj+ 1/2 = w j+\ — Wj, and so the complete artificial viscosity behaves like a second 
difference. In order to improve this situation we need to use a MUSCL or SLIP approach 
to get the first difference to higher order accuracy. Since the fluxes in the CUSP scheme 
have different slopes on the two sides of the sonic line, the SLIP construction, which relies 
on successive differences, can cause a loss of smoothness at the sonic line. This difficulty can 
be avoided by using the MUSCL approach to obtain the states corresponding to higher order 
accuracy. To impose monotonicity we can apply the limiter of the SLIP scheme. In particular, 
we can replace Eq. (12) by 


vc t 


d j+l/2 - -7r( w R 


W L ) + | [/(wr) - f(w L )] , 


where wr,Wl are given by Eq. (8). This procedure was followed throughout the numerical 
examples shown below. Application of Eq. (8) to the — (p pu pH) T variables still allows 
total enthalpy to be preserved in the higher order scheme. 
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Analysis of CUSP Scheme 

The eigenvalues of (3Arl + ucl in Eq. (14) are ft i c, /x 2 c, and /X3C. Using the simplifications of 
Eq. (22) the eigenvalues are: 

Mi = \M\, 

\M\ if \M\ < \ 

fj-2 = < a + (3 if i < M < 1 
k |M + 1| if \M\ > 1, 

\M\ if \M\ < \ 

P3 = i a — (3 if \< M < 1 
_ |M — 1| if \M\ > 1 

We note that for \M\ < ^ all three eigenvalues of the artificial viscosity are equal, and so we 
have the equivalent of a scalar viscosity. The scalar viscosity now scales with \M\c rather than 
(\M\ + l)c for the JST scalar viscosity. This is more similar to the case of preconditioning 
where ail the eigenvalues are approximately \M\c for low speed flow. Hence, we expect that 
the CUSP dissipation should work properly for very low Mach numbers provided the central 
flux terms are augmented by a suitable preconditioning matrix. 

In the subsonic range where (3 = 0, all of the versions of the CUSP scheme do not satisfy 
Eqs. (17) and (18), which are necessary for shock capturing. That is, the cell-face Mach 
numbers in the shock structure have to be larger than about 0.5 in order to avoid post-shock 
oscillations. The motivation to design uc = |u| for (3 = 0 has already been discussed previously 
and is not repeated here. However, the choice of the function for /?, as given in Eq. (19), is 
not necessarily optimal. For example, (3 = max(0,(u + \X~)/(u — A - )) would allow shock 
capturing for Mach numbers down to about 1/3, but the subsonic dissipation would be twice 
as large, uc = 2|it| for 3 = 0. Nevertheless, our own experience gained from a number of 
numerical applications suggests that there is no need for further modifications of (3. 

It is rather difficult to compare the effect of the parameter of the JST and the combined 

CUSP/ SLIP schemes since these schemes also include eigenvalue information which is not the 
same. To isolate the effect we consider a low Mach number flow with preconditioning (see 
Ref. [25] for details). Now both switches are based on the convective eigenvalue u. A typical 
value for the JST scheme is k- 4 ) = However, for a aspect ratio of one the Martinelli scaling 
[13] adds another factor of two. The parameter a = 0 in the preconditioning adds an additional 
factor of approximately 2.6. Hence, the effective constant multiplying the fourth difference is 
about which is somewhat smaller than the \ used with the original CUSP scheme with the 
SLIP formulation. For transonic flows it is more difficult to compare the levels of dissipation. 
However, it seems that the original SLIP scheme yields too high a viscosity level and so 
should be reduced to less than |. Numerical computations demonstrate the improved accuracy 
(though slower convergence) for standard transonic turbulent flows when is reduced. 

Numerical Results 

In this section we assess the accuracy and shock capturing capabilities of the HCUSP scheme. 
Comparisons are made between the HCUSP and MATD schemes. The commonly used scalar 
dissipation scheme is also included in some of the comparisons. In so doing one can clearly see 
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Figure 2a: Inviscid pressure distributions computed with MATD scheme (NACA 0012 airfoil, 
M — 0.80, a = 1.25°). 


the superiority of the high-resolution HCUSP and MATD schemes on even relatively coarse 
meshes (i.e., 16 cells in the boundary layer of a viscous flow). The flow problems considered 
in the evaluation of these numerical diffusion schemes include the following: 1) Inviscid flow 
over airfoils, 2) laminar flow over a flat plate, 3) turbulent flow over an airfoil, 4) inviscid 
and viscous hypersonic flow over a 2-D wedge. The computational effort and convergence 
behavior in computing these solutions is given. In all cases a five-stage Runge-Kutta scheme 
in conjunction with the convergence acceleration techniques of local time stepping, implicit 
residual smoothing, and multigrid was used. 

The first case is similar to the application published in Ref. [7]. Results for inviscid transonic 
flow over the NACA 0012 airfoil are shown in Fig. 1. The free-stream Mach number is 0.8 and 
the angle of attack is 1.25°. An O-topology mesh of 160 x 32 cells was used and is partially 
depicted in Fig. 1. The shock waves on the upper and lower surfaces are captured with a single 
interior point. The differences between the accurate dissipation coefficients of Eqs. (19) and 
(20) on the one hand and the simplified coefficients of Eq. (22) on the other are small. Careful 
examination reveals that the exact formulation is somewhat more accurate on the upper surface 
where the shock is stronger. These inviscid solutions were obtained with a 4-level multigrid 
and sawtooth-type cycle. They converged at a rate of about 0.90 per multigrid cycle. 

In Figs. 2 and 3 the HCUSP and MATD schemes are compared for this case. Solutions 
were computed on three successively finer C-topology meshes. The coarsest mesh contained 
192 x 32 cells, with 160 cells on the airfoil, and for each sequential mesh the number of cells 
in each coordinate direction was doubled. The principal differences between the solutions 
occur at the shock waves. Since the MATD scheme uses a pressure switch for all the flow 
equations, it cannot capture a shock with a single interior point. It requires three interior 
points. Nevertheless, the resolution of the stronger upper surface shock is nearly the same 
for both the MATD and HCUSP schemes on the 384 x 64 and 768 x 128 meshes. With the 
MATD formulation there is some smearing on the 192 x 32 mesh. As clearly evident in Fig. 3 
the HCUSP scheme allows a sharp definition of the weak lower surface shock and the Zierep 
singularity that immediately follows. The aerodynamic coefficients calculated with the MATD 
scheme on the 384 x 64 mesh essentially agree with those for the finest grid. The lift and 
drag coefficients determined with the HCUSP scheme on the corresponding meshes are not the 
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Figure 2b: Inviscid pressure distributions computed with HCUSP scheme (NACA 0012 airfoil, 
M = 0.80, a = 1.25°). 


same, with the finest grid values approaching those obtained with the MATD scheme. These 
results probably reflect the higher level of background dissipation with the HCUSP scheme. 

As an initial evaluation of the dissipation schemes for viscous flows we consider low-speed 
(Mqo = 0.15) flow over a flat plate at zero incidence. For this flow the Reynolds number per 
unit length is 10 5 . The computational domain is a rectangle. With respect to the leading edge 
of the plate, the domain extends two plate lengths upstream and one plate length downstream. 
The upper boundary is four plate lengths above the plate. Solutions were computed on the 
same domain and grids used in Ref. [23]. Starting with the finest mesh, coarser meshes were 
determined by successively eliminating every other mesh line. The finest grid consists of 
512 x 128 cells, with 384 cells on the plate. In the direction y normal to the plate the grid is 
spaced uniformly in the boundary-layer coordinate 7] (rj = y/Re x ), where x is the coordinate 
parallel to the surface, and Re x is the Reynolds number based on distance from the leading 
edge of the plate). Thus, there is constant resolution of the boundary layer at each location 
along the plate. Outside the boundary layer the grid is stretched exponentially. In order to 
resolve the region in the vicinity of the stagnation point, the grid is clustered at the leading 
edge of the plate. At the surface of the plate no-slip and adiabatic boundary conditions are 
enforced. Along the boundary upstream of the leading edge, a symmetry condition is applied. 
Characteristic type boundary conditions are used at the upstream, downstream, and upper 
boundaries. 

A comparison of the velocity profile at X/L = 0.82 computed with the the scalar, matrix, 
and HCUSP dissipation forms is displayed in Fig. 4. Even with just eight points in the boundary 
layer (64 x 16 grid) the MATD and HCUSP schemes nearly replicate the Blasius solution. As 
demonstrated in Ref. [1] scalar dissipation can produce serious contamination. With the scalar 
dissipation, more than 32 points are required in the boundary layer to obtain a grid converged 
solution. For the MATD and HCUSP schemes the variation of the errors (relative to the Blasius 
solution) in the calculated skin friction, displacement thickness, and momentum thickness are 
shown in Figs. 5a and 5b. The standard definitions [18] are used for these boundary- layer 
quantities. The errors in all the boundary-layer parameters are quite similar for the high- 
resolution schemes. This is not surprising since both schemes have a scaling factor that vanishes 
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Figure 3a: MATD solutions near lower surface shock (NACA 0012 airfoil, M = 0.80, a 
1.25°). 



Figure 3b: HCUSP solutions near lower surface shock (NACA 0012 airfoil, M = 0.80, a 
1.25°). 
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4a: Tangential and transverse velocity profiles, X/L = 0.82, 64 x 16 grid. 



U/U^ ( V/U J*SQRT (Re x ) 

4b: Tangential and transverse velocity profiles, X/L = 0.82, 128 x 32 grid. 



U/U„ (V/UJ*SQRT(Re x ) 

4c: Tangential and transverse velocity profiles, X/L — 0.82, 256 x 64 grid. 

Figure 4: Boundary-layer profiles on flat plate with M = 0.15 and Re — 10 5 . 
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Figure 5a: Comparison of MATD scheme results with Blasius solution ( M = 0.15 and Re 
10 5 ). 



Figure 5b: Comparison of HCUSP scheme results with Blasius solution (M = 0.15 and Re 

10 5 ). 
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Figure 6a: Comparison of pressure distributions with SCALAR, MATD, and HCUSP schemes 
on 160 x 32 grid (RAE 2822 airfoil, M = 0.73, a = 2.79°, Re = 6.5 x 10 6 ). 
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Figure 6b: Comparison of skin-friction distributions with SCALAR, MATD, and HCUSP 
schemes on 160 x 32 grid (RAE 2822 airfoil, M = 0.73, a = 2.79°, Re = 6.5 x 10 6 ). 
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as the surface is approached. 

We next consider transonic flow over the RAE 2822 airfoil. The free-stream Mach number 
is 0.73, the angle of attack is 2.79°, and the Reynolds number based on the airfoil chord is 

6.5 x 10 6 . Transition from laminar to turbulent flow is fixed at the 3% chord location. The 
C-type grids used are as follows: (1) 160 x 32 with 128 cells on the airfoil, (2) 320 x 64 with 
256 cells on the airfoil, and (3) 640 x 128 with 512 cells on the airfoil. The outer boundary is 
located 20 chords from the airfoil. The normal spacing at the surface of the 640 x 128 mesh is 

7.5 x 10 -6 chords. At the leading and trailing edges of the airfoil the mesh is clustered, giving 
tangential spacings of 1.17 x 10 -3 and 1.86 x 10 -3 chords, respectively. These spacings are 
roughly doubled with each mesh coarsening. To determine the effect of further mesh refinement 
a calculation was performed with the MATD scheme on a 1280 x 256 grid. 

In Figs. 6-8 the pressure ( C p ) and surface skin-friction (Cf) distributions computed with 
the different dissipation schemes for the range of meshes described are shown, along with the 
experimental data [2]. As in the inviscid cases the primary differences in the solutions occur 
at the shock wave. On the coarsest grid (160 x 32) both the SCALAR and HCUSP schemes 
produce a solution with the shock too far upstream. This is an unexpected result for the 
HCUSP scheme. The acceleration of the flow upstream of the shock is underpredicted relative 
to the finest grid. In Ref. [21] the adverse effect of a smooth limiter on the accuracy of the 
solution in the vicinity of flow transition, and thus on the acceleration of the flow upstream 
of the shock, is demonstrated. Therefore, such a result with the HCUSP scheme could be 
a consequence of the smooth limiter being used. With both the SCALAR and the MATD 
schemes a nonphysical increase in the skin-friction solution on the upper surface appears at 
the trailing edge of the airfoil. This behavior does not occur in the solution obtained with the 
HCUSP scheme. The computed aerodynamic coefficients, including the pressure and friction 
contributions to the total drag, are given in Table 1. On each mesh the lift and drag coefficients 
corresponding to the solution obtained with the MATD scheme exhibit the closest agreement 
with the 1280 x 256 grid values. There are only small discrepancies in the coefficients associated 
with the MATD and HCUSP schemes on the 640 x 128 grid. 

Convergence behavior for the HCUSP and MATD schemes is displayed in Fig. 9. Five levels 
of multigrid were used for this case and either 50 or 70 cycles were executed on two coarser 
meshes in order to obtain an initial solution. On the 320 x 64 grid the average rate of reduction 
of the residual for both schemes is about 0.92. Figure 9 also shows the effect of the modification 
given by Eq. (20) to uc in the HCUSP scheme. The convergence is improved by using the 2-D 
formulation for the dissipation coefficient uc. Note that convergence with £ = 0 was possible 
for this transonic case but not for the hypersonic case presented below. 

The fourth case is hypersonic 2-D flow over a blunt wedge. Figure 10 displays the solutions 
obtained for viscous and inviscid flow using identical meshes of 64x48 cells. Physical diffusion is 
so large that the shock profile is significantly smeared in the viscous result. For inviscid flow, on 
the other hand, we obtain perfect capturing with a single interior point in the shock structure by 
using the formulation of Eqs. (19) and (20). Detailed comparisons of the hypersonic wedge flow 
solutions yielded by the CUSP scheme and AUSM have been presented in Ref. [14]. The shock 
capturing capabilities of both schemes are essentially equal. A comparison of shock profiles 
for the exact and the simplified coefficients is given in Fig. 11. We have chosen the first-order 
scheme in order to address the pure shock capturing capability of the CUSP scheme without 
interference from the limiter. The simplified dissipation coefficients of Eq. (22) produce strong 
oscillations at the shock, even though there is substantial physical diffusion present. Hence it is 
concluded that an accurate implementation of dissipation coefficients is a must for hypersonic 
flows with strong shocks. 
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Comparisons of computation times indicate that the CUSP scheme needs about 40% more 
computer time than the basic scalar dissipation. The MATD scheme only requires about 
15% additional time. MATD requires less CPU time primarily because it needs only a single 
evaluation of the limiter function. Due to lower inherent dissipation, computations with the 
CUSP formulation converge somewhat slower for transonic flows than those with simple scalar 
dissipation. The major advantage of the CUSP approach is that it is more accurate and 
more robust than scalar viscosity. Our numerical tests indicate that the accuracy of the 
CUSP scheme is intermediate between scalar and matrix dissipation for transonic flows. For 
hypersonic flows it seems to be more robust than the matrix viscosity even though it is not 
TVD. 

Since the CUSP scheme is implemented through artificial dissipative terms, it does not 
have to be applied at each stage of the Runge-Kutta method. In particular, the diffusive fluxes 
can be evaluated only at the first, third, and fifth stages of a five-stage method, just as is 
typically done for the scalar numerical dissipation. 


Concluding Remarks 

The CUSP scheme has been studied and analyzed. A detail comparison has been made between 
the CUSP, MATD, and scalar dissipation schemes. For transonic inviscid flows the CUSP 
scheme allows better resolution of shock waves, since they are captured with one interior 
point. However, the aerodynamic quantities such as lift and drag obtained with the CUSP 
scheme are not as accurate on coarser meshes (i.e., 320 x 64 cells or less) as those calculated 
with the MATD scheme. Both the CUSP and MATD formulations give high accuracy in 
the computation of high Re number flat-plate flow. For transonic viscous flows and coarser 
meshes the accuracy in aerodynamic coefficients is better with the MATD scheme than with 
the CUSP scheme. This loss in accuracy with the CUSP scheme on coarser grids appears to be 
a consequence of limiter activation (i.e., reduction to first order). Convergence of the HCUSP 
scheme has been improved by introducing the aspect-ratio scaling factor of Eq. (20) . 

With our present choice of HCUSP dissipation coefficients it has been shown that the 
resolution of strong shock waves occurring in hypersonic flows is possible whereas the simplified 
coefficients that were published previously failed. 

The CUSP scheme requires roughly 40% more computer time than the scalar scheme, while 
the MATD scheme needs about 15% more time. Using only pressure in the switches for MATD 
rather than using a different measure for each equation reduces the robustness of the algorithm 
but requires less computer time. Convergence behavior with the CUSP and MATD schemes 
is similar. For hypersonic flows the CUSP scheme seems to be more robust than the MATD 
scheme. 
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Figure 7a: Comparison of pressure distributions with SCALAR, MATD, and HCUSP schemes 
on 320 x 64 grid (RAE 2822 airfoil, M = 0.73, a = 2.79°, Re = 6.5 x 10 6 ). 



Figure 7b: Comparison of skin-friction distributions with SCALAR, MATD, and HCUSP 
schemes on 320 x 64 grid (RAE 2822 airfoil, M — 0.73, a = 2.79°, Re = 6.5 x 10 6 ). 
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Figure 8a: Comparison of pressure distributions with SCALAR, MATD, and HCUSP schemes 
on 640 x 128 grid (RAE 2822 airfoil, M = 0.73, a = 2.79°, Re = 6.5 x 10 6 ). 



Figure 8b: Comparison of skin-friction distributions with SCALAR, MATD, and HCUSP 
schemes on 640 x 128 grid (RAE 2822 airfoil, M = 0.73, a = 2.79°, Re = 6.5 x 10 6 ). 
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Dissipation 

Scheme 

Grid 

c L 

C D 

Cd p 

Cdj 

SCALAR 

160 x 32 

0.8172 

0.01728 

0.01275 

0.004532 


320 x 64 

0.8331 

0.01743 

0.01194 

0.005487 


640 x 128 

0.8532 

0.01782 

0.01225 

0.005574 

MATD 

160 x 32 

0.8304 

0.01818 

0.01251 

0.005662 


320 x 64 

0.8538 

0.01808 

0.01250 

0.005571 


640 x 128 

0.8597 

0.01799 

0.01246 

0.005535 

| 

1280 x 256 

0.8611 

0.01800 

0.01246 

0.005544 

HCUSP 

160 x 32 

0.7987 

0.01926 

0.01367 

0.005594 


320 x 64 

0.8493 

0.01831 

0.01263 

0.005679 


640 x 128 

0.8592 

0.01803 

0.01245 

0.005585 


Table 1: Lift and drag coefficients for turbulent flow over RAE 2822 airfoil (RAE 2822, M x , = 0.73, 
a = 2.79° , Re c = 6.5 x 10 6 ). 
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9b: Effect of modified vc on HCUSP scheme. 


Figure 9: Convergence histories for turbulent flow computations (RAE 2822 airfoil, = 0.73, 
a = 2.79°, Re c = 6.5 x 10 6 ). 
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center line flow 



Figure 11: Influence of HCUSP dissipation coefficients on hypersonic flow over 2-D wedge. 
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